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1.  Introduction 


The  Tonge-Ramesh  material  model  is  a  micromechanics  based  constitutive  model 
for  the  high  rate  failure  of  quasi-brittle  materials  such  as  armor  ceramics.  It  incorpo¬ 
rates  micromechanics  based  damage,  granular  flow  of  the  damaged  material,  lattice 
plasticity  (volume  preserving),  and  equation  of  state  coupling.  The  key  physical  as¬ 
pects  of  the  model  are  described  in  (Tonge  2014,  eh.  4).  This  document  serves  to 
detail  the  user  input  parameters  and  additional  functionality  that  is  not  discussed 
in  the  paper.  This  document  is  not  intended  to  discuss  the  physical  implications  of 
the  input  parameters  or  discuss  the  physical  reasoning  used  to  develop  the  material 
model. 

2.  Material  Model  Inputs 

The  following  list  describes  the  model  inputs.  All  inputs  must  be  provided,  and  the 
current  model  does  not  do  any  error  checking  to  ensure  that  reasonable  values  are 
provided  for  the  user  inputs.  Units  are  described  using  the  following  abbreviations: 

•  P :  Stress — Unless  otherwise  noted,  scalar  measures  of  stress  and  strain  used 
in  this  model  are  the  magnitude  of  the  deviatoric  part  of  the  tensor  and  the 
magnitude  of  the  isotropic  part  of  the  tensor.  For  a  generic  tensor  A,  r a  = 
\/Aiev  :  Aiev  and  za  =  tr(A)/\/3.  The  same  measures  of  deviatoric  and 
isotropic  magnitude  are  used  for  both  stress  and  strain.  Constants  from  other 
models  that  use  traditional  measures  like  the  Von-Mises  equivalent  stress  and 
the  pressure  will  need  to  be  rescaled. 

•  L:  Length 

•  r:  Time 

•  T :  Temperature 

•  M:  Mass 

•  E:  Energy 

1)  useDamage:  Unitless — Flag  to  activate  damage  calculation  module;  0 — do 
not  perform  damage  calculation;  1 — do  damage  calculation. 
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2)  usePlasticity:  Unitless — Flag  to  activate  J2  plasticity;  0 — do  not  per¬ 
form  calculation;  1 — perform  calculation. 

3)  useGranularPlasticity:  Unitless — Flag  to  activate  pressure  depen¬ 
dent  granular  flow  calculations;  0 — do  not  perform  calculation;  1 — perform 
calculation. 

4)  useOldStress:  Unitless — Flag  that  controls  sublooping  in  damage  growth 
calculation;  1 — do  not  do  sublooping,  instead  use  the  stress  from  the  begin¬ 
ning  of  the  time  step. 

5)  artif  icialViscosity:  Unitless — Flag  to  activate  linear  and  quadratic 
bulk  viscosity  calculation;  0 — do  not  perform  calculation;  1 — perform  calcu¬ 
lation. 

6)  art  if  icialViscousHeat  ing:  Unitless — Flag  to  convert  work  done  by 
artificial  viscous  pressure  to  temperature  rise;  0 — do  not  perform  calculation; 
1 — perform  calculation. 

7)  BulkModulus:  P — Isothermal  bulk  modulus  at  zero  pressure.  Must  be  con¬ 
sistent  with  the  equation  of  state  (K0  =  pC%). 

8)  ShearModulus:  P — Shear  modulus  (G0  in  Eq.  2). 

9)  rho:  ML~ 3 — Reference  material  density  (p0). 

10)  FlowStress:  P — Flow  stress  for  J2  plasticity  (r0  in  Eq.  12). 

1 1)  HardeningModulus:  P — Isotropic  linear  hardening  modulus  for  J2  plas¬ 
ticity  ( H  in  Eq.  12). 

12)  InitialPlasticStrain:  Unitless — Initial  plastic  strain  for  J2  plastic¬ 
ity  calculation.  This  is  a  repetitive  input  and  its  meaning  could  change  in 
future  versions. 

13)  J2RelaxationTime:  r — Timescale  for  visco-plastic  relaxation  of  ^plas¬ 
ticity. 

14)  NumCrackFamilies:  Integer — Number  of  flaw  bins  or  families  used  to 
discretize  the  flaw  size  distribution  (Ayms  in  Eq.  10). 

15)  MeanFlawSize:  L — Mean  flaw  size  used  only  for  flaw-distribution  types: 
0  (delta)  and  1  (Gaussian). 
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16)  FlawDensity:  L  3 — Average  number  of  flaws  per  unit  volume. 

17)  StdFlawSize:  L — Standard  deviation  of  the  flaw  size  distribution.  Used 
only  in  distribution  type  1  (Gaussian). 

18)  FlawDistType:  Integer — Distribution  type  for  the  probability  distribution 
of  flaw  sizes;  0 — Delta;  1 — Gaussian;  2 — Pareto.  The  majority  of  the  model 
development  has  been  done  using  a  Pareto  distribution. 

19)  MinFlawSize:  L — Minimum  flaw  size  for  the  flaw  size  distribution. 

20)  MaxFlawSize:  L — Maximum  flaw  size  for  the  flaw  size  distribution. 

21)  FlawDistExponent:  Unitless — Exponent  for  the  Pareto  flaw  size  distri¬ 
bution. 

22)  RandomizeFlawDist:  Unitless — Flag  to  apply  flaw  size  distribution  ran¬ 
domization  during  the  first  step  of  the  simulation;  1 — apply  the  randomiza¬ 
tion;  0 — do  not  apply  the  randomization. 

23)  RandomSeed:  Integer — One  part  of  the  random  seed  used  to  initialize  the 
random  number  generator  for  setting  the  initial  flaw  size  distribution. 

24)  RandomizeMethod:  Integer — Set  the  method  used  to  randomize  the  flaw 
size  distribution.  Choices  are  1-7;  7  is  the  recommended  choice,  while  3  will 
cause  each  successively  finer  flaw  bin  to  have  twice  the  flaw  density  until  all 
of  the  flaw  bins  are  filled.  When  using  Method  3,  the  smallest  flaw  may  be 
smaller  than  MinFlawSize  depending  on  the  size  of  the  element  volume. 
Method  3  may  be  an  approach  to  use  fewer  bins  to  cover  the  flaw  distribu¬ 
tion.  See  (Tonge  2014,  app.  B)  for  a  complete  description  of  the  different 
randomization  methods. 

25)  BinBias:  Unitless — Scale  the  discretization  of  the  flaw  distribution  to  better 
capture  small  flaws.  Only  used  when  RandomizeMethod  =  6  or  7. 

26)  KIc:  PL0,5 — Fracture  toughness  for  the  micromechanics-based  damage  cal¬ 
culation  (Kjc  in  Eq.  11). 

27)  FlawFriction:  Unitless — Coefficient  of  friction  for  clean  self  contact. 
Used  in  the  micromechanics  calculation. 
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28)  FlawAngle:  Radians — Angle  between  the  maximum  compressive  stress 
and  the  flaw  normal.  See  (Tonge  2014,  fig.  4.3)  for  an  illustration. 

29)  FlawGrowthExponent:  Unitless — Exponent  for  the  micromechanics  crack 
growth  rate  (yc  in  Eq.  11). 

30)  FlawGrowthAlpha:  Unitless — Denominator  for  reference  crack  growth 
rate  in  micromechanics  calculation  (ac  in  Eq.  11). 

31)  Crit  icalDamage:  Unitless — Damage  level  to  activate  granular  plasticity. 

32)  MaxDamage:  Unitless — Maximum  allowable  damage  level.  Once  damage 
reaches  this  value,  the  damage  calculation  is  no  longer  performed. 

33)  MicroMechPlaneStrain:  Unitless — Flag  to  select  plane  strain  (1)  or 
plane  stress  (0)  for  the  micromechanics  based  damage  calculation  (recom¬ 
mend  1). 

34)  MaxDamagelnc:  Unitless — Maximum  increment  in  damage  if  damage  based 
sublooping  is  enabled. 

35)  UseDamageTimestep:  Unitless — flag  to  allow  the  damage  growth  rate  to 
determine  the  global  time-step.  This  is  a  legacy  parameter  that  is  no  longer 
used.  It  should  be  set  to  0. 

36)  dt_increaseFactor:  Unitless — Increase  factor  for  the  time-step  when 
using  the  damage  growth  rate  to  compute  the  time-step  size.  This  is  a  legacy 
parameter  that  is  no  longer  used.  It  should  be  set  to  100. 

37)  IncInitDamage:  Unitless — Flag  to  include  initial  flaw  size  in  the  damage 
calculation;  0 — do  not  include  the  initial  flaw  size  in  the  calculation  of  the 
material  point  damage;  1 — compute  damage  from  both  the  initial  flaw  size 
and  the  wing  crack  size. 

38)  DoFlawInteract  ion:  Unitless — Flag  to  activate  the  self  consistent  inter¬ 
action  calculation.  This  calculation  can  be  a  significant  portion  of  the  damage 
growth  calculation  time.  This  may  be  a  possible  place  to  trade  accuracy  for 
speed;  0 — skip  the  calculation  and  use  the  element  stress  for  the  crack  growth 
calculation;  1 — perform  the  self  consistent  calculation. 

39)  GPTimeConst:  r — Relaxation  time  for  granular  flow. 
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40)  JLoc:  Unitless — Inelastic  volume  change  ratio  to  activate  the  localized  flag. 

Must  be  greater  than  1. 

41)  GP  Granular  SI  ope:  Unitless  if  GPYieldSurf  Type  =  1,  P  if  GPYieldSurf  Type 

=  2 — Slope  of  the  granular  flow  surface.  Parameter  A  in  Eq.  13  and  Eq.  14. 

This  input  has  units  of  P  when  used  with  GPYieldSurf  Type=2. 

42)  GP  Cohesion:  P — Intersection  of  the  granular  flow  surface  with  the  hydro¬ 
static  stress  axis  (parameter  B  in  Eq.  13  and  Eq.  14). 

43)  GPYieldSurf  Type:  Unitless — Switch  between  a  cone  (1)  and  parabola 
(2)  type  yield  surface  using  the  2  surface  representation.  The  parabola  (2)  is 
not  recommended. 

44)  GPPc:  P — Full  consolidation  pressure  when  using  the  2  surface  granular 
flow/pore  compaction  model. 

45)  GP  Jr  e  f :  Unitless — Reference  volume  change  ratio  for  the  2  surface  granular 
flow/pore  compaction  model. 

46)  GPP  ref:  P — Reference  pressure  for  the  2  surface  granular  flow/pore  com¬ 
paction  model. 

47)  GPSurf  aceType:  Unitless — Flag  that  controls  the  use  of  the  single  surface 
granular  flow  model  (1)  or  the  2  surface  (any  value  other  than  1)  model  used 
in  the  original  material  model  development.  The  single  surface  model  is  under 
active  development  and  may  change  in  future  releases.  It  is  not  recommended 
for  production  runs. 

48)  AbsToll:  Unitless — Absolute  tolerance  for  single  surface  granular  flow  re¬ 
turn  calculation  (recommend  1  x  1CT14). 

49)  RelToll:  Unitless — Relative  tolerance  for  single  surface  granular  flow  re¬ 
turn  calculation  (recommend  1  x  10“8). 

50)  Maxlter:  Unitless — Maximum  iterations  for  single  surface  granular  flow 
return  calculation  (recommend  20). 

51)  MaxLevels:  Unitless — Maximum  number  of  recursive  refinement  levels 
for  granular  flow  calculation.  (This  functionality  can  dramatically  increase 
the  model  cost  and  does  not  significantly  improve  the  granular  flow  update 
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accuracy.  We  recommend  that  this  parameter  be  set  to  1  to  disable  automatic 
refinement.) 

52)  GFMSmO:  Unitless — Parameter  /i0  in  Eq.  22. 

53)  GFMSml:  Unitless — Parameter  /p  in  Eq.  22. 

54)  GFMSm2:  Unitless — Parameter  /i2  in  Eq.  22. 

55)  GFMSpO:  Unitless  (less  than  0) — Parameter  p0  in  Eq-  23.  Saturation  com¬ 
paction  mean  stress  at  large  values  of  evp.  This  is  the  mean  stress  required  to 
begin  compressing  out  porosity  for  the  fully  distended  material.  This  value 
should  be  almost  0. 

56)  GFMSpl:  P — Parameter  p1  in  Eq.  23. 

57)  GFMSp2:  Unitless — Parameter^  in  Eq.  23. 

58)  GFMSp3:  Unitless — Parameter^  in  Eq.  23. 

59)  GFMSp4:  Unitless  (0  <  <  1) — Parameter^  in  Eq.  24. 

60)  GFMSal:  P — Parameter  a±  in  Eq.  20.  Hydrostatic  intercept  for  linear  portion 

of  Ff. 

61)  GFMSa2:  P_1 — Parameter  a2  in  Eq.  20.  Nonlinear  portion  of  Ff  inside  the 
exponent. 

62)  GFMSa3:  P — Parameter  a.>  in  Eq.  20.  Exponential  scale  for  Ff. 

63)  GFMSBeta:  Unitless  (0  <  (3  <  1.0) — Scale  value  to  introduce  nonassocia- 
tive  behavior  into  the  bulking  response  of  the  material. 

64)  GFMSPs  i :  Unitless — Ratio  of  the  triaxial  tensile  strength  to  the  triaxial  com¬ 
pressive  strength  (Brannon  et  al.  2009). 

65)  GFMSJ3Type:  Unitless — Select  the  form  for  the  Lode  angle  dependence 
(Brannon  et  al.  2009);  0 — Drucker-Prager;  1 — Gudehus;  2 — William- Warnke. 

66)  ArtViscl:  Unitless — Linear  artificial  bulk  viscosity  coefficient. 

67)  ArtVisc2:  Unitless — Quadratic  artificial  bulk  viscosity  coefficient. 
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68)  MGCO:  L/t — Bulk  sound  speed  at  reference  conditions  for  Mie-Griineisen 
equation  of  state. 

69)  MGGammaO:  Unitless — Griineisen  coefficient  for  equation  of  state. 

70)  MGS1:  Unitless — Slope  of  Us  —  Up  curve. 

71)  MGS  2:  Unitless — Parameter  for  equation  of  state. 

72)  MGS  3:  Unitless — Parameter  for  equation  of  state. 

73)  MGCv:  fiM"* 1 2 3/!”1 — Specific  heat  per  unit  mass. 

74)  MGTheta:  T — Reference  temperature. 

75)  JMin:  Unitless — Unused. 

76)  MGNPts:  Unitless — Unused. 


Not  all  constants  are  used  in  all  computations.  Constants  10-13  are  only  used 
if  usePlast icity  =  1.  Constants  15-38  are  only  used  if  useDamage  =  1. 
Within  these  constants,  15  and  17  are  only  used  if  FlawDistType  =  0  or  1. 
Constants  21  and  25  only  apply  when  FlawDistType  =  2.  Constants  39-65 
only  apply  when  useGranularPlasticity  =  1.  Constants  41-46  apply  when 
GPSurf  aceType  =  0.  Constants  48-65  apply  when  GPSurf  aceType  =  1. 

3.  Material  Model  State  Variables 

This  material  model  computes  the  evolution  of  the  total  Cauchy  stress  (the  devia- 
toric  and  isotropic  parts)  following  the  Abaqus  UMAT  conventions.  In  addition  to 
updating  the  Cauchy  stress,  the  model  tracks  the  following  internal-state  variables. 
They  are  stored  in  an  array,  but  are  given  names  if  the  host  code  allows  (this  is  not 
a  standard  capability  in  an  Abaqus  UMAT). 

1)  Iel  ( [I ei  =  |tr (bei))  :  Unitless — Must  be  greater  or  equal  to  1.0. 

2)  damage  (D):  Unitless — Must  be  positive  and  is  no  greater  than  MaxDamage. 

3)  JGP  (JGp  ~  exp  ^ tr(dGp)dr^):  Unitless — Inelastic  volume  change  ratio 
associated  with  granular  flow. 
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4)  GP_s train  (7 GP  =  f*  \\ddGP  || dr):  Unitless — Accumulated  deviatoric  in¬ 
elastic  granular  flow. 

5)  GP_energy  (f*  Jcr  :  dGPdr ):  EL~ 3 — Energy  dissipated  by  granular  flow 
per  unit  reference  volume. 

6)  plasticStrain  (yp  =  f*  \ \dp\ \dr):  Unitless — Equivalent  plastic  strain 
for  the  .U  type  plastic  flow. 

7)  plasticEnergy  (f*  jp\ |  Jcrdev| \dr):  EL~3 — Energy  dissipated  through  J2 
plasticity  per  unit  reference  volume. 

8)  artViscPres:  P — Artificial  viscous  pressure  that  should  be  added  to  the 
stress  tensor. 

9)  volHeating:  TV-1 — Rate  of  change  of  the  temperature. 

10)  localized:  Unitless — Flag  that  is  set  when  localization  conditions  are 
met:  1 — Damage  is  greater  than  or  equal  to  Crit  icalDamage;  2 — Granular 
flow  localization  conditions  met  (JGP  >  JLoc);  3 — Conditions  for  1  and  2 
were  both  met;  4 — Localization  caused  by  J2  flow  (the  computed  strength 
due  to  hardening  is  0);  5 — Conditions  1  and  4  were  met;  6 — Conditions  2 
and  4  were  met;  7 — All  3  localization  conditions  were  met  (1,  2,  and  4). 

11)  epsVGP  ( evGP  =  f*  tr(dGp)dr):  Unitless — Inelastic  volume  strain.  Only 
used  when  GPYieldSurf  Type=l. 

12)  gamGP  (7 Gp  =  /„  WdGP  ||  dr):  Unitless — Accumulated  Inelastic  deviatoric 
strain.  Only  used  when  GPYieldSurf  Type=l. 

13)  epsVGP_qs:  Unitless — Same  as  epsVGP,  but  for  the  rate-independent  so¬ 
lution.  This  is  only  used  when  the  single  surface  model  (GPYieldSurf  Type=l) 
is  used  with  rate -dependent  granular  flow  (GPTimeConst>  0).  The  visco¬ 
plasticity  algorithm  (Brannon  2007)  tracks  both  the  actual  material  state  and 
the  state  that  represents  the  rate-independent  solution. 

14)  gamGP_qs:  Unitless — Same  as  gamGP,  but  for  the  rate-independent  solu¬ 
tion.  This  is  only  used  when  the  single  surface  model  (GPYieldSurf  Type=l) 
is  used  with  rate -dependent  granular  flow  (GPTimeConst>  0).  The  visco¬ 
plasticity  algorithm  (Brannon  2007)  tracks  both  the  actual  material  state  and 
the  state  that  represents  the  rate-independent  solution. 
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15)  si g_qs_l  1 :  P — Stress  associated  with  the  rate-independent  solution  to  the 
granular-flow  evolution.  This  is  only  used  when  the  single  surface  model 
(GPYieldSurf Type=l)  is  used  with  rate-dependent  granular  flow 
(GPTimeConst>  0).  The  visco-plasticity  algorithm  (Brannon  2007)  tracks 
both  the  actual  material  state  and  the  state  that  represents  the  rate-independent 
solution. 

16)  si g_qs_2 2 :  P — Stress  associated  with  the  rate-independent  solution  to  the 
granular-flow  evolution.  This  is  only  used  when  the  single  surface  model 
(GPYieldSurf Type=l)  is  used  with  rate-dependent  granular  flow 
(GPTimeConst>  0).  The  visco-plasticity  algorithm  (Brannon  2007)  tracks 
both  the  actual  material  state  and  the  state  that  represents  the  rate-independent 
solution. 

17)  s ig_qs_3 3:  P — Stress  associated  with  the  rate-independent  solution  to  the 
granular-flow  evolution.  This  is  only  used  when  the  single  surface  model 
(GPYieldSurf Type=l)  is  used  with  rate-dependent  granular  flow 
(GPTimeConst>  0).  The  visco-plasticity  algorithm  (Brannon  2007)  tracks 
both  the  actual  material  state  and  the  state  that  represents  the  rate-independent 
solution. 

18)  si g_qs_2 3 :  P — Stress  associated  with  the  rate-independent  solution  to  the 
granular-flow  evolution.  This  is  only  used  when  the  single  surface  model 
(GPYieldSurf Type=l)  is  used  with  rate-dependent  granular  flow 
(GPTimeConst>  0).  The  visco-plasticity  algorithm  (Brannon  2007)  tracks 
both  the  actual  material  state  and  the  state  that  represents  the  rate-independent 
solution. 

19)  s  ig_qs_l  3:  P — Stress  associated  with  the  rate-independent  solution  to  the 
granular-flow  evolution.  This  is  only  used  when  the  single  surface  model 
(GPYieldSurf Type=l)  is  used  with  rate-dependent  granular  flow 
(GPTimeConst>  0).  The  visco-plasticity  algorithm  (Brannon  2007)  tracks 
both  the  actual  material  state  and  the  state  that  represents  the  rate-independent 
solution. 

20)  s  i g_qs_l  2 :  P — Stress  associated  with  the  rate-independent  solution  to  the 
granular-flow  evolution.  This  is  only  used  when  the  single  surface  model 
(GPYieldSurf Type=l)  is  used  with  rate-dependent  granular  flow 
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(GPTimeConst>  0).  The  visco-plasticity  algorithm  (Brannon  2007)  tracks 
both  the  actual  material  state  and  the  state  that  represents  the  rate-independent 
solution. 

21)  f lawNumber_k:  L~ 3 — Flaw  number  density  in  bin  number  k  (c Ok). 

22)  starterFlawSize_k:  L — Initial  flaw  size  in  bin  number  k  (sfe). 

23)  wingLength_k:  L — Flaw  growth  in  bin  number  k  (4). 


Unused  state  variables  are  set  to  0.  The  material  model  also  updates  the  temperature 
of  the  material  in  the  TEMP  function  argument  for  the  Abaqus  UMAT  interface. 

4.  Key  Physical  Equations 
4.1  Elastic  Response 

The  model  assumes  a  decoupled  representation  of  the  Kirchhoff  stress  tensor: 

T  =  Tdev  -  psJeI.  (1) 


The  deviatoric  stress  Tdev  is  a  linear  function  of  the  deviatoric  part  of  the  volume 

—  _ 2/3  T 

preserving  elastic  deformation  as  measured  by  be  =  Je  FeF e  : 


T  dev 


G 


(2) 


Here  G  is  the  damaged  shear  modulus  defined  as 


G(D)  = 


2D 

15" 


(3  Zr  +  2  Zn 


(3) 


The  parameters  Zn,  Zr,  and  Zc  are  functions  of  the  elastic  moduli  and  relate  to  the 
compliance  of  an  individual  crack. 


The  volumetric  response  is  determined  by  a  Mie-Gruneisen  equation  of  state: 

+  PoFc)  [ec(Je)  +  cv(9  —  0O)]  (4) 

if  JP  <  1 


Ps(Je,»)=IljDPH(je 

An 


1  -  y  (1  -  Je) 


P0Cn2(l-Je 


pH(Je)  =  )  (l-5l(l-Je)-S2(l-Je)2-53(l-Je)3)2 

PqCq(1  —  Je)  otherwise 
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(5) 


The  scaling  factor  K/K0  softens  the  bulk  modulus  in  response  to  damage  growth. 
The  damaged  bulk  modulus  is  defined  as 


K{D)  =  (Kq1  +  D  (Zn  +  4 Zc))  1 .  (6) 


Here  K0  is  the  undamaged  bulk  modulus  (— dps/dJe\D=0).  The  cold  energy  ec  is 
the  solution  to 


dec 
— -  + 
dJe 


eo) 


Ph_ 

Po 


r 

2 Je 


(1  -  Je) 


(7) 


Here  eo  =  cv9o  and  the  initial  condition  is  ec(l)  =  0.  At  finite  temperature  (90  >  0) 
the  solution  for  ec  is  not  strictly  positive.  Since  the  reported  strain  energy  SSE  is 
ec  added  to  the  deviatoric  component  ( G(D )  (tr(Bel)  —  3)  /2),  the  reported  strain 
energy  may  not  be  positive  for  all  deformations. 


The  artificial  viscosity  formulation  is  the  linear  quadratic  viscosity  commonly  used 
in  many  hydrocodes: 


Pvisc 


(AiC0\tr(de)\dx  +  A2lr(de)2dx2)  if  tr(de)  <  0 


otherwise 


(8) 


The  heating  rate  associated  with  the  work  done  by  the  artificial  viscosity  is 


9 vi  sc 


'  J Pviscd’(de 
PoCy 


(9) 


4.2  Micromechanics  of  Damage  (useDamage  =  1) 

The  model  uses  a  micromechanics-based  damage  model  where  damage  is  defined 
as 

-Nbins 

D  =  Uk  (sfc  +  /fc)3  .  (10) 

k= 1 

In  the  summation,  the  bin  number  k  loops  over  the  Abins  that  are  used  to  discretize 
the  local  flaw  size  distribution;  u>k  is  the  number  density  of  flaws  per  unit  volume 
that  are  represented  by  the  flaw  family  k\  the  initial  flaw  size  is  s^,  which  has  grown 
an  additional  length  lk  due  to  the  applied  loading  history. 

When  DoFlawInteract ion  =  1,  the  model  uses  a  self-consistent  Eshelby  in- 
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elusion  solution  to  compute  the  elevated  stress  in  the  neighborhood  of  a  flaw  due 
to  the  surrounding  flaws.  This  solution  involves  the  solution  to  a  linear  system  of  8 
equations  and  can  be  disabled  (at  the  expense  of  not  including  the  interactions)  by 
setting  DoFlawInteract ion  =  0.  After  computing  the  effective  stress  around 
a  flaw  (creq),  the  stress-intensity  factor  for  crack  growth  is  computed  based  on  a 
wing-cracking  mechanism  (Tonge  2014,  eqs.  4.57  and  4.58).  It  depends  on  the 
crack-face  coefficient  of  friction  (FlawFrict  ion),  applied  load  (creq),  the  flaw 
angle  (FlawAngle),  and  the  representative  flaw  size  .sg.  From  that  stress-intensity 
factor  (A/),  we  compute  the  crack-growth  rate  using 


cv  (  K  f  I<lc  Vc 
ac  U/  —  0.5A ic ) 


(ID 


The  increased  crack  length  is  used  to  update  the  damage  parameter. 


4.3  Traditional  J2  Plasticity  (usePlasticity  =  1) 

The  J2  volume  preserving  plasticity  implementation  uses  a  multiplicitive  decompo¬ 
sition  of  the  deformation  gradient  into  elastic  and  plastic  portions.  It  is  integrated  in 
the  spatial  configuration  using  an  objective  integration  algorithm  (Simo  and  Hughes 
2000,  ch.  9).  Activation  of  this  component  may  not  be  justified  for  all  ceramic  sys¬ 
tems.  It  is  an  elastic -plastic  material  model  with  linear  strain  hardening.  The  yield 
surface  is  defined  by 

/(t)  =  11'7-devil  -  (“qP)  (HeP  +  T°)  •  (12) 

The  factor  22  corrects  for  the  effect  of  damage  growth  where  the  effective  yield 
stress  decreases  with  damage  in  the  same  way  that  the  shear  modulus  decreases 
with  damage  (the  elastic  strain  required  to  cause  yield  is  independent  of  damage). 
The  hardening  modulus  may  be  either  positive  or  negative  (softening).  When  the 
material  softens  to  zero  strength,  the  localized  flag  is  set  to  4.  The  coupling 
between  this  plasticity  module  and  damage  is  not  heavily  tested. 


4.4  Granular  Flow  (useGranularPlasticity  =  1) 

Two  separate  granular-flow  models  are  implemented  in  this  material  model  version. 
The  original  model  was  discussed  in  (Tonge  2014,  ch.  4)  and  uses  a  2-surface  ap¬ 
proach  to  describe  the  distortion  and  compaction  mechanisms  for  granular-flow.  A 
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newer  single-surface  model  has  been  added  and  is  based  on  the  tear-drop-shaped 
single-surface  model  (Fossum  and  Brannon  2006;  Brannon  et  al.  2009). 

4.4.1  Two-Surface  Flow  Compaction  Model  (GPSurfaceType  =  0) 

This  2-surface  model  computes  the  granular-flow  and  porosity  increase  first  fol¬ 
lowed  by  the  porosity  compaction  calculation.  The  compaction  calculation  does 
not  depend  on  the  deviatoric  stress. 


4.4.1. 1  Granular  Flow  and  Porosity  Generation  Flow  behavior  is  associa¬ 
tive  to  the  yield  surface  and  there  is  no  hardening  and  this  flow  formulation  is 
written  in  terms  of  the  Kirchhoff  stress  r.  There  are  2  choices  for  the  shape  of 
the  granular-flow-yield  surface  determined  by  the  input  GPYieldSurf  aceType. 
Surface  1  is  defined  by 

-a).  (13) 

Yield  Surface  2  changes  the  units  of  A  to  pressure  units  and  is  defined  by 

-  b)  •  (14) 

A  linear  viscosity  model  is  activated  by  setting  GPTimeConst  to  a  positive  value. 
The  linear  viscosity  model  follows  a  Deviant-Louis-type  viscoplasticity  model  that 
is  unconditionally  stable  and  only  tracks  the  current  state  of  the  material  (Simo  and 
Hughes  2000,  sec.  2.7.4). 

The  plastic  flow  direction  is  m  =  (n  +  (AI) / y/l  +  (2A2.  Here  n  =  rdev/|  |rdev|  |, 
1  =  1/ y/3,  A  is  the  tangent  to  the  yield  surface  in  r  —  z  space,  and  (  is  a  parameter 
that  is  iteratively  decreased  from  1  to  ensure  positive  plastic  dissipation. 


/(t)  =  rdev  :  rdev  +  A 


tr(r 


f{r)  =  \J t dev  :  rdev  +  A 


tr(r 

VI 


4.4.1 .2  Pore  Compaction  This  module  is  turned  on  when  granular  plasticity  is 
activated  using  the  2-surface  model  (GPSurf aceType=0).  This  is  an  additional 
yield  surface  that  depends  on  only  the  hydrostatic  pressure  (p  =  —  l/3tr(cr)). 
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It  is  defined  by 


ft(p,jGF,j)  =  i 


Pc-Po 

(. JGP 


P 


P-Po  (  tGP  _  tGP\ 


P<Po 

Po<P  <  Pc 


[  JGP  ~  1  P<PC- 

(15) 

This  is  a  simple  porosity  model  with  a  quadratic  crush  behavior  for  pressures  be¬ 
yond  P0  and  an  exponential  compaction  behavior  for  lower  pressures.  The  plastic 
flow  resulting  from  pore  compaction  is  purely  isotropic. 


4.4.2  Single  Surface  Flow  Compaction  Model  (GPSurfaceType  =  1) 

This  granular-flow  model  is  based  on  the  Kayenta  yield  and  flow  surface  (Brannon 
et  al.  2009).  The  yield  surface  is  written  in  terms  of  the  Lode  invariants.  For  a 
Cauchy  stress  tensor  cr  with  a  deviatoric  part  s  =  cr  +  pi,  the  Lode  invariants  are 


r  —  y/s  :  s 

(16) 

z  =  — =tr  (cr) 

(17) 

.  J3  /  3  \  2/3 

3^tU  ' 

(18) 

The  yield  function  (Fossum  and  Brannon  2006)  is  written  as 


f(r,z,0,%,e;)  =  r(0)2ir2  -  F,(z,  %) \F,(z,lr)\ Fc(z,e ;) 
F/i vj  —  f"i  -  <i3exp  (a2VSz)  -  HoVlz) 


Fc(z,e°)  = 


1  -  if  P3z  < 

(K(e*p)-X(e%)) 


K 


otherwise 


M7p)  =  A* 1  +  (a*o  -  Pi)  exp  (-/i27 p) 
x(eP  =  Pi  (po  +  ^  (ln(p3  +p2evp)  -  |ln(p3  +p26p)|) 
n(evp)  =  p4X  (el). 


(19) 

(20) 

(21) 

(22) 

(23) 

(24) 


The  plastic  flow  direction  is 


M  —  a(N  +  f3N  iso). 


(25) 
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Here  N dev  and  N iso  are  the  deviatoric  and  isotropic  parts  of  yield  surface  normal, 
/ 3  is  an  input  parameter  (  fi  <  1),  and  a  is  a  normalization  parameter  to  ensure  that 
M  :  M  =  1  (Brannon  and  Leelavanichkul  2010,  eq.  21).  This  granular-flow  model 
uses  the  multistage  return  algorithm  (Brannon  and  Leelavanichkul  2010). 

5.  Yield  Surface  Illustrations 

For  illustrative  purposes,  the  micromechanics-based,  damage-growth  model  can  be 
viewed  as  providing  an  evolving  damage-growth  surface  that  depends  on  the  ex¬ 
treme  principal  stresses.  Such  a  surface  is  illustrated  in  Fig.  1.  For  a  delta  distribu¬ 
tion  of  cracks,  crack  growth  begins  when  the  stresses  exceed  the  blue  line.  Under 
compressive  loading,  the  stress  required  to  sustain  damage  growth  increases  until 
the  stability  limit  (green  line),  then  decreases  with  increasing  damage  until  a  user- 
specified  maximum  allowable  damage  level  is  reached.  Figure  1  shows  the  stress 
history  for  a  material  point  loaded  at  a  constant  strain  rate  under  uniaxial  stress 
conditions.  The  material  point  stress  is  plotted  as  a  function  of  time.  The  peak  in 
the  blue  curve  represents  the  stability  limit  for  this  loading  condition  and  the  transi¬ 
tion  to  the  dotted  green  line  represents  the  activation  of  granular  flow  when  damage 
reaches  the  user-defined  threshold  of  0.125. 


Most  Compressive  Principal  Stress  (GPa)  Damage 

a)  b) 

Fig.  1  a)  shows  the  compressive  stress  states  beyond  which  damage  growth  occurs,  with  the  ar¬ 
row  indicating  the  path  for  uniaxial  compression;  b)  plots  uniaxial  compressive  stress  required 
for  damage  growth  as  a  function  of  damage  level. 


Once  granular  flow  is  activated,  there  is  a  competition  between  porosity  growth  and 
pore  compaction.  The  pressure  required  to  initiate  pore  compaction  as  a  function  of 
the  distension  (a  measure  of  porosity)  is  shown  in  Fig.  2.  The  2  different  choices  for 
GPSurfaceType  use  different  representations  of  the  crush  curve,  but  either  can 
be  fit  to  experimental  data  if  available. 
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Crush  curves 


Fig.  2  Illustration  of  the  crush  curves  showing  pressure  required  to  initiate  pore  collapse  for 
pure  hydrostatic  loading  as  a  function  of  distension  for  the  2  different  granular-flow  surface 
types 


The  single- surface  representation  for  the  granular-flow-yield  surface  is  enabled  by 
setting  GPSurf  aceType=l.  This  model  provides  a  smooth  representation  of  the 
increase  in  deviatoric  strength  with  pressure  until  pore  compaction  begins  to  dom¬ 
inate,  then  provides  a  gradual  reduction  in  strength  with  increasing  pressure  up  to 
the  crush  pressure.  The  meridional  profile  of  the  yield  surface  is  shown  in  Fig.  3. 
This  model  includes  the  option  to  specify  Lode  angle  dependence,  which  reduces 
the  strength  of  the  material  as  the  loading  conditions  deviate  from  triaxial  compres¬ 
sion.  For  different  ratios  of  triaxial  tensile  strength  to  triaxial  compressive  strength, 
the  octahedral  profile  is  shown  in  Fig.  4. 
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Fig.  3  Meridional  profile  of  the  yield  surface  for  granular  flow  when  using 

GPSurf aceType=l 


Fig.  4  Octahedral  profile  of  the  yield  surface  for  granular  flow  when  using 

GPSurf  aceType=l  and  GFMS  J3Type=2  for  different  values  of  the 
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6.  Single  Material  Point  Driver 


The  source  package  for  the  material  model  comes  with  a  single  element  driver  for 
the  material  model  called  runUMATTR_MixedBC.  This  small  program  drives  the 
material  model  using  a  user-defined  combination  of  stress-  and  strain-rate  boundary 
conditions  as  a  function  of  time.  This  driver  is  called  using  this  command  line: 


. /runUMATTR_MixedBC  <Boundary  File>  <Parameter  File>  <History  File> 


The  Boundary  File  is  a  text  file  containing  white-space-separated  numbers. 
The  first  line  contains  the  time- step  size  and  the  representative  edge  length  followed 
by  6  flags  (0  or  1)  where  0  denotes  an  rate  of  deformation  (symmetric  part  of  the 
velocity  gradient)  and  1  denotes  an  applied  Cauchy  stress-boundary  condition  for 
each  of  the  symmetric  tensor  components  in  the  order  XX  YY  ZZ  YZ  XZ  XY.  The 
subsequent  lines  contain  a  time  followed  by  the  boundary  conditions.  The  applied 
boundary  conditions  are  treated  as  piecewise  constant;  they  are  not  interpolated 
between  time  points  in  the  table.  The  simulation  ends  at  the  last  time  in  the  table. 

The  Parameter  File  is  a  text  file  that  contains  all  of  the  material  parameters 
listed  in  Section  2  one  parameter  per  line  with  the  tag  followed  by  white  space  fol¬ 
lowed  by  the  value  for  that  parameter.  The  History  File  is  an  output  file  with 
column  headers  on  the  first  line  then  fixed-width  columns  of  all  history  variables  at 
every  simulation  time  step  in  the  subsequent  lines. 


18 


7.  Model  Test  Suite 


The  granular-flow  algorithm  is  tested  against  3  analytic  solutions  (Brannon  and 
Leelavanichkul  2010).  The  first  of  3  analytic  solutions  tests  an  associative  Drucker- 
Prager  material  subjected  to  plastic  loading  that  causes  a  rotation  of  the  principal 
directions  of  the  stress  tensor  while  maintaining  constant,  principal  stress  values. 
The  results  from  this  verification  test  are  shown  in  Fig.  5. 


Brannon  Ex  1 
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Fig.  5  Depiction  of  the  loading  path  in  a  stationary  3-dimensional  manifold  for  stress  space 

The  second  problem  tests  a  J2  Von-Mises  material  under  plastic  loading  that  is  not 
coaxial  to  the  stress  deviator.  A  graphical  comparison  of  the  simulation  results  with 
the  analytic  solution  is  provided  in  Fig.  6. 
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Fig.  6  Comparison  between  the  simulated  stress  components  and  the  analytic  solution  under 
nonproportional  loading  for  a  Von-Mises  material 


The  third  test  problem  is  a  5-stage  loading  path  for  an  associative  Drucker-Prager 
material  expanded  from  the  published  version  (Brannon  and  Leelavanichkul  2010) 
to  include  a  leg  that  traverses  around  the  hydrostatic  tensile  vertex.  The  expanded 
solution  for  associative  flow  was  provided  by  M.  Hommel.  The  first  stage  is  hydro¬ 
static  compression.  The  next  leg  of  the  problem  is  triaxial  extension  loading  (the  2 
equal  principal  values  are  greater  than  the  third)  such  that  continued  loading  after 
yield  results  in  a  stationary  stress.  The  loading  direction  is  then  changed  so  that  the 
stress  crosses  the  yield  surface  and  approaches  the  yield  surface  in  the  normal  direc¬ 
tion.  Upon  yield,  the  stress  travels  along  the  yield  surface  to  higher  deviatoric  and 
hydrostatic  stresses.  The  next  leg  is  plastic  loading  to  the  hydrostatic  tensile  vertex. 
After  reaching  the  hydrostatic  tensile  vertex,  plastic  loading  takes  the  stress  back  to 
the  location  where  yield  occurred  during  the  second  leg  of  the  test.  The  results  from 
the  simulation  are  shown  in  Figs.  7  and  8.  The  results  do  not  exactly  match  the  an¬ 
alytic  solution;  however,  component  tests  that  isolated  the  granular-flow  portion  of 
the  model  demonstrated  significantly  better  agreement  when  the  bulk  modulus  and 
shear  modulus  were  exactly  constant  and  the  additive  decomposition  of  the  strain 
increment  was  consistent  with  the  total  strain  measure.  The  use  of  a  Mie-Griineisen 
equation  of  state,  the  coupling  between  the  bulk  modulus  and  the  distension,  and  the 
geometric  nonlinearity  introduced  by  the  use  of  bel  to  define  the  deviatoric  stress  all 
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contribute  to  the  differences  between  the  analytic  solution  and  the  simulation  results 
shown  in  Figs.  7  and  8. 


Fig.  7  Comparison  between  the  simulated  stress  components  and  the  analytic  solution  for  the 
5-leg  loading  path 


r 


Fig.  8  Illustration  of  the  simulated  and  input  stress  path  in  a  signed  deviatoric  stress  versus 
hydrostatic  stress  space 


Verification  of  the  material  model’s  implementation  in  the  host  code  was  performed 
by  testing  simple  shear  at  a  shear  rate  of  10,000  1/s  in  all  6  directions  (12,  13,  21,  23, 
31,  23).  The  stress  history  from  a  representative  simple  shear  simulation  is  shown 
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in  Fig.  9.  Shear  in  the  other  6  directions  produced  similar  results.  All  other  history 
variables  show  similar  agreement  among  the  different  host  codes. 


Fig.  9  Stress  history  for  representative  single-element  simulations  showing  multiple  host  codes 
producing  nearly  identical  results  for  simple  shear  to  a  shear  of  1 


8.  Model  Interface  Notes 


The  model  uses  an  Abaqus  UMAT  interface,  which  has  a  function  prototype  of 


void  PTR_umat_stressUpdate ( 

double  STRESS [6],  double  STATEV [ ] ,  double  DDSDDE [ 6 ] [ 6 ] , 
double  *SSE,  double  *SPD,  double  *SCD,  double  *RPL, 
double  DDSDDT [ 6] ,  double  DRPLDE [ 6 ] ,  double  *DRPLDT, 

const  double  STRAN[6],  const  double  DSTRAN[6],  const  double  TIME [2], 

const  double  *DTIME,  double  *TEMP,  double  *DTEMP, 

const  double  *PREDEF,  const  double  *DPRED,  const  double  *CMNAME, 

const  int  *NDI,  const  int  *NSHR,  const  int  *NTENS,  const  int  *NSTATV, 

const  double  PROPS [],  const  int  *NPROPS,  const  double  COORDS [3], 

const  double  DR0T[3] [3],  double  *PNEWDT,  const  double  *CELENT, 

const  double  DFGRD0[3] [3],  const  double  DFGRD1[3] [3], 

const  int  *NOEL,  const  int  *NPT,  const  int  *LAYER, 

const  int  *KSPT,  const  int  *KSTEP,  const  int  *KINC) ; 


22 


The  following  exception  types  may  be  thrown  by  the  constitutive  model  if  some 
invalid  conditions  are  encountered  (like  negative  plastic  work): 

•  runt ime_error 

•  out_of_range — should  not  happen.  Indicates  an  incorrect  matrix  index. 

•  domain_error — Matrix  math  error,  either  inversion  of  a  singular  matrix 
or  division  by  zero. 

Currently  all  implementations  allow  thrown  exceptions  to  propagate  and  cause  the 
host  code  to  crash. 

Total  internal  energy  can  be  computed  from  SSE+MGCv*  (TEMP-MGTheta). 

DROT  is  the  incremental  rotation  that  is  consistent  with  the  stress-integration  al¬ 
gorithm.  STRAN  is  an  approximation  to  the  logarithmic  strain,  but  in  this  model 
only  the  3  normal  terms  are  used  to  compute  the  volume-change  ratio.  To  compute 
the  total  strain,  one  can  integrate  the  strain  rate  using  the  same  objective  algorithm 
used  to  integrate  the  stress.  Total  plastic  work  per  unit  of  mass  is  SPD.  This  is  the 
same  as  (plast icEnergy+GP_energy)/rho.  Axisymmetric  geometry  is  not 
supported  because  it  is  inconsistent  with  the  failure  modes  in  brittle  materials.  The 
longitudinal  sound  speed  is  sqrt  (DDSDDE  [  0  ]  [  0  ]  /rho_cur)  . 

The  function  PTR__umat_stressUpdate_ALE3D  tracks  and  updates  the  tem¬ 
perature  in  the  SCD  variable  place  holder.  This  is  a  work-around  for  host  codes  that 
do  not  accurately  track  the  temperature  being  fed  into  the  UMAT. 
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